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Abstract 

We present approximative solutions to the neutrino evolution equa- 
tion calculated by different methods. In a two neutrino framework, 
using the physical parameters which gives the main effects to neutrino 
oscillations from v e to another flavors for L ~ 3000Km and E ~ lGeV 
, the results for the transition probability calculated by using series 
solutions, by to take the neutrino evolution operator as a product of 
ordered partial operators and by numerical methods, for a linearly 
and sinusoidally varying matter density are compared. The extension 
to an arbitrary density profile is discussed and the evolution operator 
as a product of partial operators in the three neutrino case is obtained. 

1 Introduction 



The conversion of one type of neutrino into another, while propagating in vac- 
uum or in matter, can be the solution to the Solar and atmospheric neutrino 
problems and to explain the results about neutrino flux from experiments 
with reactors and accelerators. In order for such vacuum or matter oscil- 
lations to occur, it is necessary that the neutrino have a nonzero mass and 
that the neutrino masses be not all degenerate. These possibilities have been 
explored in several models like that ones where a term of mass for neutrinos 
is constructed by the introduction of new Higgs fields and the mass is gen- 
erated by spontaneous breakdown of symmetry, models with fermion sector 



1 



extended and more complex models. For a brief review about these mod- 
els, see [1]. In this context flavour eigenstate of neutrinos can be a mixing 
of mass eigenstates and oscillations can take place. Independent of specific 
models, for a neutrino initially created as a flavour eigenstate, the amplitude 
of probability for conversion to another flavor and the survival amplitude of 
probability, after to travel a distance L, in vacuum or in matter [2], are given 
by the solutions of the evolution equation. If CP conservation is assumed the 
evolution equation depends on five parameters, three mixing angles and two 
differences of squared masses. There is exact solution to the evolution equa- 
tion for neutrinos propagating in vacuum or in a constant matter density in 
a two neutrino framework [3]. The three neutrino case in a constant matter 
density has been analyzed in the references [4], [5], [6]. The exact expression 
to the evolution operator for a three neutrino system in a constant matter 
density in terms of the mass squared differences and the vacuum mixing an- 
gles has been firstly derived in the reference [5], the effective mixing angles 
in matter have been also derived. In the reference [6] an alternative expres- 
sion has been derived. For varying matter density solutions can be found 
in the adiabatic approximation [7 ] and in the nonadiabatic case near the 
resonances [8]. These solutions can give a estimation about the flux of neu- 
trinos in the cases which the neutrino source is big and there is a continuum 
spectrum of energy, the integral on the energy or on the point of neutrino 
production cancels contributions from the interference terms. For situations 
where the dimension of the source is shorter than the oscillation lengths, for 
determined neutrino energy, numerical solutions for the evolution equation 
can give informations about mixing angles and differences of squared masses. 
This is the case for neutrinos from reactors and accelerators and atmospheric 
neutrinos travelling across the earth at energy E m lGeV . The evolution 
equation for a neutrino system propagating at variable matter density has 
recently been studied . In the reference [9], in a two neutrino framework, 
the strategy is to expand the neutrino oscillation probability and the effec- 
tive potential in series. In the reference [10] the strategy of Elisio Lisi and 
Daniele Montanino is to take the evolution operator given by the ordered 
product of partial evolution operators from the shells along of the neutrino 
trajectory. In many of the three flavour descriptions of neutrino oscillations 
the so called 'one mass scale dominance' 5m 2 = Ara^ << Am^ ~ Am| 2 
has been adopted. In these cases, for particular situations, the problem re- 
duce to the two neutrinos one. This is the case of the reference [11] where 
the simple approach Ui 3 ~ 1 has been used and in the ref. [12] where the 
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influence of the mixing angle 6i 3 has been tacked into account for Sub-GeV 
atmospheric neutrinos and the accuracy of the approximation is determined 
by the relation tg29\ 3 ~ 29\ 3 . In these cases evolution equations have been 
derived and the solutions in a varying matter case needs to be found with the 
help of some approximative method. In ref. [13], in a two neutrino scenario, 
approximate solution to the evolution equation has been presented for high 
energetic neutrinos, E > 10 GeV for neutrinos propagating at the matter of 
the Earth. In ref. [14] approximate solution to the neutrino equation in a 
two neutrino framework with corrections due to the three neutrino mixing 
has been presented for low energetic neutrinos , the validity of the solution 
depends on the parameter e = 2 ■ The main problem of to take the evo- 
lution operator as a ordered product of partial operators in the three neutrino 
case is that the expression to the evolution operator in this case is not simple. 
In reference [15] this strategy has been adopted but the approach Am^ = 
reduce in some sense the complications to that ones from the two neutrino 
case. In reference [16] the general case has been analyzed using the results 
from [6]. In this work, in section II, in order to obtain numerical solutions by 
methods like Runge Kuta, the evolution equation of the neutrino system at 
variable electron density is written as a system of two disconnected second 
order linear ordinary differential equations. The conversion probability is 
computed for neutrinos travelling at linear matter density. In order to com- 
pare our results with another ones the parameters of the system are those 
used in reference [9] . In section III the evolution operator of the neutrino 
system at variable electron density is computed as the product of partial 
operators. The conversion probability is computed for the case analyzed in 
section II and we apply the methods to compute the conversion probability 
for neutrinos travelling in a sinusoidally varying matter density. In section 
IV we apply the strategy of to take the evolution operator as a product of 
ordered partial operators to the three neutrino case using the results from 
reference [5]. The conversion probabilities are computed for the constant and 
the linearly varying matter density cases Finally, in section V, we summarize 
our results and give our conclusions. 
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2 Neutrinos in matter at variable density 

In a two neutrino framework the evolution equation for a neutrino system 
propagating at matter is given by 



where 



"/(.'■)=( ?S I (2) 



H(x) = A30-3 + A\G\ and we have defined A 1 = ^A sin 26* and A 3 = 
■^jjj(A(x) — A cos 29) , (T3 and o\ are Pauli matrices and in these expressions, 
A(x) = 2^/2G F N e (x)E , A = wij — m\ denote the mass squared difference 
between the two mass eigenstates in vacuum, 9 is the mixing angle and N e (x) 
and E are the electron density and neutrino energy respectively. Introducing 
the new fields $(x) by the relation 

Uf (x) = e-^ dXAWa3 ^(x), (3) 

X 

using the auxiliary matrices o± = ai± ™ 2 , the definition f3(x) = J d\A(X), 
and with the help of the Backer-Hausdorff relation for exponential operators 

00 „ 

e~ F Ge F = V ^f; G = G; G n = [G n _i, F] (4) 

n=0 

the evolution equation, in terms of the field &(x) = ( tf l ) , reads 



i^(x) = A 1 [e l2 ^a + + e- lW{x) a_]<$>(x) (5) 

or 

i^i(z) = A ie ^Mx) (6) 

z^0 2 (x)=A 1 e-^)0 1 (x) 

solving this equations for the fields </>i(x) one obtain the following linear 
system of second order differential equations [17] 

p ; - 22^3(3;)-^ + ^(rr)0i(rr) = (7) 
4 



For neutrinos produced as u e for example, the initial conditions to the fields 
4>i(x) becomes 

MO) = 1; 0i(O) = 0; 2 (O) = 0; 0' 2 (O) = -iij (8) 

and the equations (7) with the initial conditions (8) can be solved by using 
some numerical method for an arbitrary matter density. The survival and 
conversion amplitudes of probability are given by the expressions 

u e (x) = e-^Mx) (9) 
u^x) = e^Mx) 

We apply this results for a linearly varying matter density N e (x) = N core (l + 
f)/4 where N core is the electron density in the earths' core and L = 3000 
km. In the numerical calculations, we have used the mixing angle 6 = 13° 
and the squared mass difference Am 2 = 2 x !0~ 3 eV 2 . These choices are 
made in order to our results can be compared with the results from [9]. The 
solutions has been found numerically. In Fig. 1 the transition probability as 
a function of distance is presented. 



3 The evolution operator 

The solution of (1) in a general way is Vf{x) = U(x,x )b , f(x ) where U(x,x ) 
is the evolution operator of the neutrino system given explicitly by 

U(x,x ) = Exp < -i J H(X)dX > (10) 

\ XO ) 

where H is the Hamiltonian of the system and the symbol Exp concern to 
cxpansional that represent a expansion in orderly multiple integrals 

fx \ x x A 

-i J H(X)d\ > = l-i-i J H(X)dX+(i) 2 J H(X)dX-i J H(X')dX'+... 

(11) 



XO 



XO XO XO 
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and that for [H(x),H(x')] = is reduced to usual exponential. This is the 
case of neutrinos travelling the vacuum and the matter of constant density. 
The strategy of Elisio Lisi and Daniele Montanino [10] is to divide the neu- 
trino trajectory in infinitesimal elements j in which the electronic density is 
given by Nj(x) = N](x) + 5Nj(x) 

Xj 

Nj(x) = / N e {x)dx (12) 

Xj Xj—\ J 

Xj-! 

and the neutrino propagation is described by the Hamiltonian H(x) = Hj + 
SHj(x). In this way, the evolution operator is taken as an product of ordered 
partial operators 

U(x,x ) =Y[U j (x j ,Xj- 1 ) (13) 

j 

Note that for a linearly varying matter density the quantities Nj(x) and 
Hj(x) correspond to N e (x) and H(x) where x = . 

The Hamiltonian Hj(x), except for the proportional term to the identity 
that just introduces in the evolution operator a constant phase factor is given 
by 

„, s If V^GpN) - kcos29 ksm26 1 £TT . N 

2l ksnTB -V2G F W J + kcos29 j +6H ^ (M) 

where k = and 6Hj(x) = ^°£™M a3 . 
The operator in the jth shell is 

Uj{xj, x^) = Exp i-i J [7fj + 5Hj(x)]dx > . (15) 

[ X 3-l ) 

With the help of the relation between exponential operators [18] 

Exp \ -i J [F(X) + G(X)]d\ > =Exp\-i J F(X)dX \ * (16) 



Exp I -i j E 1 (X)G(X)E 2 (X)dX 
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where 



E'i(A) = Exp |-i j F(u)du 
E 2 (\) = Exp j-i y F(u)du| 



and by note that if., is constant, the evolution operator Uj(xj, Xj_i) becomes 

f 7 - 

(17)' 

and this expression permit us to take the evolution operator in orders of 
perturbation. At the first order Uj(xj,Xj-i) results 

Xj-l 

where the first term corresponds to the operator for a neutrino travelling 
at a constant matter density. The expression to the zero order term can be 
obtained easily but, in order to prepare the method to the three neutrino 
case, when this expression is more complicated, we summarize the steps for 
to do this. The way is to calculate the eigenvalues A12 to the Hamiltonian Hj 
, Ai, 2 = ^ where w = VA 2 + A 2 - 2AA cos 29 and to obtain the expression 
to the evolution operator in the particle basis. The result is 



ttP , 



U j (a^Xj-i) = j Q e - ia /AE j ( 19 ) 

Then, the evolution operator in the flavor basis is given by 



rr I cos 9 sin 9 -^p cos 6* — sin 9 1 

Uj{Xj,Xj-{)= \ ~ n)Uj\ . ~ ~ ) (20) 

— sin 9 cos 9 J sine* cos 9 ' 



where 

~ A sin 29 
sin 29 = 

w 
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and 9 denotes the effective mixing angle in matter. With the definitions 

k m = 2]jj > c j = cos[k m (xj — Xj_i)/2] , Sj = sin [/c m (xj — Xj-i)/2], Uj(xj, Xj-i) 
becomes [10] 

yj , s c, + i-Sj cos 29 —is,-, sin 29 , 

Uj{Xj, Xj-i) — < J J ~ J ~ >. (21) 

[ — isj sin 2$ Cj — isj cos 2# 1 

The first order term demand a more hard work and we will explicit only the 
final expression 

6U(x j ,x j - 1 ) = - t -Bm2el C > Sin J° C 3 cos29-iS 3 
y 3 3 ' 2 \ C jC os29 + iSj -C^ sin 20 ' 



where 



Sj = V2Gp J dxSNj(x) sin k m (x — x) 



Xj-! 



J 

Cj = \p2G F J dx5Nj(x) cosk m (x — x) 



Xj-! 



and x = (xj — Xj-i)/2. 

The partial evolution operator when taken in first approach order is not 
automatically unitary, being necessary to normalize it. That unitarity is 
guaranteed by the redefinition 



U?(x) = Ei+^M (23) 

where 



I ^4 + [S] + Cj] sin 2 29 



B 2 

By the utilization of the evolution operator as a product of ordered 
infinitesimal operators (13) we calculate the transition probability as a func- 
tion of distance in zeroth-order approximation for the same case presented 
in section II, a linearly varying matter density. The results are presented in 
Fig.l. We calculate the same quantity by utilization of the method used in 
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ref. [9] , the result is presented in Fig. 1 also for comparison. For neutri- 
nos travelling in matter with an arbitrary density profiles we must to take 
into account that not necessarily Nj(x) = N e (x). As an illustrative exam- 
ple we apply the method for a sinusoidally varying matter density N e (x) = 
N core (l + smnx). With the approximation sm £^ x ~ 1 where Ax is the length 
of the shells in the neutrino trajectory the relation Nj(x) = N e (x) can be 
used. In Fig. 2 the transition probability as a function of distance for neutri- 
nos at energy 500MeV and the physical parameters used in previous sections 
is presented. The oscillation length is ^ = 4000/cm. The same quantity 
calculated by using Runge-Kuta method is presented also for comparison. 



4 The three neutrino case 

In the three neutrino case de evolution equation becomes 

% Tx Vf[x) = 2^E j™ 2 ^" 1 + 2 } v * x ) ( 24 ) 
where M 2 is the 3x3 diagonal matrix of mass 

Vf (x) = (u e (x), u,(x), u T (x) f ; U = j*M e i*M e iu>A a (25) 

and A is the matrix which the only different from zero element is An = A(x). 
Defining the quantities A = m\ — m\ ; E = m| + m\ ; A x = £ — 2m\ 
and A = E — A c2uj, for a constant matter density, in terms of the field 
<&(#) = e~^ A? 'Vf(x) the evolution equation becomes 

d . . ~ . . 
i—$(x) = H $(x). 

(jjJb 



where 



or 



(A 31 + A 32 )c 2 A31 + A32 Ac2ujc 2 (j) 2A 



where 5H is the relevant term for transitions which elements are given by ^ 



^ = 2 + 3 2~ + T 
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A 3 i + A 32 A c2u A _ A s2u c<p 
h 22 - g + - 3 ; h 12 

(A 31 + A 32 )s 2 A31 + A32 Ac2a>s 2 A 
h3 ~ 2 + 3 2 3 

A s2cu s0 (A 3 i + A 32 + A c2uj)s2(j) 

«23 ^ ' ^ 13 ~~ 4 • 

The eigenvalues of the hamiltonian if are the effective quantities ^| in pres- 
ence of matter and are related to de roots \ of the hamiltonian 5H by the 

relation _ 

mf £ + ml + A 

% _ 2 4- X 

2E 6E + l - 

With the help of the quantities Q and R 

f A 2 A? A 2 AAc 2 0c2u AAi(c 2 (p - 2s 2 0) 

Q = -{- + ii + T 2^ + ^^ — - 

-1 \A\ 9A X A 2 3AA 2 (c 2 0-2s 2 0) 3A 2 A!(c 2 - 2s 2 0) 

il = — — \ — : h 



27 t 4 4 4 2 

-1 / 9A 2 Ac 2 0c2cu 9AA 2 (c 2 - 2s 2 0) VAA^^u 3 
27 \ 2 4 2 

cos a = 



27 



the roots A,- results 



-1 -Q a 

Ax = ——\ cos — 

E V 3 3 



1 -Q a 1 r — . a 

A 2 = — \ / cos v —Q sin — 

2E\ 3 3 2£ V V 3 



. 1 /-<2 a 1 r -— . a 

A3 = — \ / cos — v —Q sin — . 

3 2E\ 3 3 2E V ^ 3 

The elements of the symmetric evolution operator T(x,o ) for the field 
<&(x) can be found in terms of these roots [5] and result 

T ^ = ^ [(Aj - h 22 )(Xj - h 33 ) - j4| p . a . l/2E 



d 



jkl 
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r 22 = Z 



3 




r 12 



E 



3 




(26) 



2.3 



E 



3 



[{Xj - h u )h 23 + hi 2 h 13 ] iXiX /2E 
djki 
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E 



[(Aj - h 22 )h 13 + h u h 23 ] e -i\ jX /2E 
djki 



where djki = (Xj — Xk)(Xj — A;), with j ^ k ^ I, and the evolution operator 
in the flavor basis equals: 



For varying matter density the strategy is the same of the expression 
(13) taking the evolution operator at zeroth order of approximation in the 
jth shell given by expression (27). For more practical numerical calculations 
is most easy to calculate the evolution of the field <&(#) when the partial 
operator is T(xj,Xj-i) and after to obtain the amplitudes Vf(x) with the 
relation between the fields Vf(x) and <&(x).In the figure [1] we present the 
transition probability P ex = P eil + P er as a function of distance in zeroth- 
order approximation in the three neutrino case (TNC) for a linearly varying 
matter density the mixing angle uj = 28.7° and the difference of squared mass 
A 2 i = 5 x lO~ 5 eV 2 are chosen from indications of Solar neutrinos [19]. In the 
figure [3] the transition probability P ex computed by using the approximate 
evolution operators (13) and the partial operators T(xj, Xj-\) given by (27) is 
presented for neutrinos at energy E = 400 MeV and the same parameters of 
the figure [1], the probabilities P 6(U and P er are presented also. In the figure[4], 
for neutrinos at energy E = 500MeV the probabilities P ex calculated by using 
expressions (13) and (27) and by using the exact solution (27) are presented 
for neutrinos travelling in a constant matter density N e (x) = 1.5 N core . The 
case uj = 7° is presented also. 



G(x,x ) = e l ^ A7 T(x,x )e-^ A7 . 



(27) 
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5 Summary and conclusions 



We have presented approximate solutions to the neutrino evolution equation 
calculated by different methods. The Fig. 1 shows the agreement between 
numerical results and that one obtained by using the strategy of Elisio Lisi 
and Daniele Montanino and a little disagreement between these results and 
that one calculated using the strategy of Mattias Blennow and Tommy Ohls- 
son for L ps 3000 km. One can see in ref [9] that the convergence of the 
transition probability calculated by using series solutions with numerical re- 
sults is better as much terms are included in these series and for different 
numbers of terms the disagreement occur at different values of L. The re- 
sults calculated by the utilization of the evolution operator as a product of 
ordered partial operators at first-order approximation do not introduce sig- 
nificative differences and are not presented. Additionally for an arbitrary 
density profile using the strategy of Lisi at zeroth-order approximation is 
necessary only to calculate the appropriate Nj(x) according to we can see in 
the figure [2] for a sinusoidally varying matter density. In contrast, by using 
the strategy of Blennow is necessary a lot of calculations for each new law 
of variation density. The differences of to calculate the approximate evolu- 
tion operators as a product of ordered partial operators in a two and three 
neutrino cases are clearly showed in the figure [3] where P ex is presented for 
neutrinos at E — 400 MeV . In the three neutrino case the calculation of 
any quantity P a p is performed directly from expression (27). In the figure [4] 
we have superposition of the probabilities calculated by using ordered partial 
operators and the exact solution to the neutrino evolution equation and the 
case of uj = 7° tend to the two neutrino case in order to demonstrate the 
consistency of the method. Meanwhile the main question from this figure 
is: For constant matter density N e (x) = 1.5 N core and the parameters used 
in this case, the relation A 3 i >> A 2 i is real but there is a disagreement 
between P ex calculated by using the two and three neutrinos scenarios. As 
the neutrino penetrate more and more in the matter, the initial little differ- 
ences will propagate along of the neutrino trajectory and the disagreement 
between these approaches will appear for great values of L. For neutrinos at 
E = 1200 MeV travelling across the Earth matter this effect will appear also 
in agreement of the results presented in the figure [1] near L ~ 3000i^m. Fi- 
nally, we have derived an approximative expression to the evolution operator 
for a system of three flavors of neutrinos propagating in matter. Differently 
of the approaches used in references [12], [13], [14] and [15] the method can 
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be used for any values of the neutrino energy, mixing angles in vacuum and 
difference of squared masses and applied to neutrinos propagating at mat- 
ter with an arbitrary density profile. The expressions to the components of 
the evolution operator in a constant matter density (27), from reference [5], 
are more compact and practical for to construct the partial operators, than 
ones from the reference [6]. Due to the form of the matrix A of interaction 
between neutrinos and matter, the extension for to introduce CP violation 
needs only a little modification in the expression (27). 
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Figure Caption 1 



13 



The neutrino oscillation probability P ex for a linearly varying matter den- 
sity as a function of the distance for E = 1.2GeV. The solid black line cor- 
responds to the numerical solution, the red line to the series solution and 
the blue line and the black dashed line are the results obtained using the 
evolution operator as a product of ordered partial operators at zeroth-order 
approximation in the two and three neutrino cases respectively. 

Figure Caption 2 

The neutrino oscillation probability P ex for a sinusoidally varying matter 
density N e {x) = N core (l + smKx) as a function of the distance. The neutrino 
energy is 500MeV and the oscillation length is ^ = 4000A;m. The solid red 
line corresponds to the numerical solution and the blue line is the results ob- 
tained using the evolution operator as a product of ordered partial operators 
at zeroth-order approximation. 

Figure Caption 3 

The neutrino oscillation probability P ex for a linearly varying matter den- 
sity as a function of the distance for E = AOOGeV. The black line and the 
red line correspond to the results obtained using the evolution operator as 
a product of ordered partial operators at zeroth-order approximation in the 
two and three neutrino cases respectively. The blue and green lines are the 
probabilities P eti and P er obtained in a three neutrino context with the angle 
tp = 0. 

Figure Caption 4 

The neutrino oscillation probability P ex for a constant matter density as a 
function of the distance for E = 500 MeV. The black line corresponds to the 
results obtained using the evolution operator as a product of ordered partial 
operators at zeroth-order approximation in the two neutrino case. The red 
line and the blue dashed line are the probabilities P ex obtained in a three 
neutrino context approximately and exactly. The green line corresponds to 
the approximative case with u = 7°. 
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